
Autor: Celeste Castro Granados $\mathbb C \hat{e}l \mathbb s$
Importamos primero la biblioteca que utilizaremos para graficar y probar nuestras funciones
using Plots #Biblioteca para graficar
Vamos a obtener la 2nda derivada de una función en el punto $x_i$ con una precisión en la que se desprecian términos en $h^5$, es decir, tenemos que realizar operaciones con dos puntos a la izquierda y dos a la derecha para una expresión simétrica (orden 4). Los puntos van a estar siempre separados por una distancia h porque se trata de una malla uniforme.
Por lo tanto, partimos de las siguientes ecuaciones (expansiones en serie de Taylor):
$$ f(x_0+h) = f(x_0) + f'(x_0)h + \frac{f''(x_0) h^2}{2} + \frac{f'''(x_0) h^3}{3!} + \frac{f^{iv} (x_0) h^4}{4!} + \theta (h^5) \quad...(1)$$$$ f(x_0-h) = f(x_0) - f'(x_0)h + \frac{f''(x_0) h^2}{2} - \frac{f'''(x_0) h^3}{3!} + \frac{f^{iv} (x_0) h^4}{4!} + \theta (h^5)\quad...(2) $$$$ f(x_0+2h) = f(x_0) + 2f'(x_0)h + 2f''(x_0) h^2 + \frac{4f'''(x_0) h^3}{3} + \frac{2f^{iv} (x_0) h^4}{3} + \theta (h^5) \quad...(3)$$$$ f(x_0-2h) = f(x_0) - 2f'(x_0)h + 2f''(x_0) h^2 - \frac{4f'''(x_0) h^3}{3} + \frac{2f^{iv} (x_0) h^4}{3} + \theta (h^5) \quad...(4) $$Sumamos las ecuaciones (1)-(2) y (3)-(4) y obtenemos lo siguiente:
$$ f(x_0+h) + f(x_0 - h) = 2f(x_0) + f''(x_0) h^2 + \frac{1}{12} f^{iv} (x_0) h^4 + \theta (h^5) \quad...(5)$$$$ f(x_0+2h) + f(x_0 - 2h) = 2f(x_0) + 4f''(x_0) h^2 + \frac{4}{3} f^{iv} (x_0) h^4 + \theta (h^5) \quad...(6)$$Así, podemos hacer 16Eq(5) - Eq(6):
$$ 16 f_{i+1} + 16 f_{i-1} - f_{i+2} - f_{i-2} = 30 f_i + 12 f''_i h^2 $$Por lo que, la expresión para la segunda derivada numérica de orden 4 en una malla uniforme es:
$$ f''(x_i) = \frac{1}{12h^2} \bigg ( 16 f_{i+1} + 16 f_{i-1} - f_{i+2} - f_{i-2} - 30 f_i \bigg ) $$Observamos que la expresión anterior no puede ser aplicada en los puntos extremos. Para calcular la derivada en los extremos necesitamos las expresiones para $x_0+h$ , $x_0+2h$, $x_0+3h$ y $x_0+4h$. Estas expresiones podemos escribirlas como el siguiente sistema de ecuaciones:
$$\begin{bmatrix} 1 & 1/2! & 1/3! & 1/4! \\ 2 & 2 & 4/3 & 2/3 \\ 3 & 9/2 & 9/2 & 27/8 \\ 4 & 8 & 32/3 & 32/3 \end{bmatrix} \begin{bmatrix} h f_i' \\ h^2 f_i'' \\ h^3 f_i ''' \\ h^4 f_i^{(iv)} \end{bmatrix} = \begin{bmatrix} f_{i+1} - f_i \\ f_{i+2}-f_i \\ f_{i+3} - f_i \\ f_{i+4}-f_i \end{bmatrix}$$$\therefore$ La solución para $f_i''$ es la solución del segundo nivel y es la siguiente:
$$ f''(x_i)= \frac{\frac{35}{12} f_i -\frac{26}{3} f_{i+1} + \frac{19}{2} f_{i+2} - \frac{14}{3} f_{i+3} + \frac{11}{12} f_{i+4} }{h^2} $$Para el primer punto se evalúa en $i=1$ y para el último en $i=end$, considerando que:
$i+1 \rightarrow i-1, i+2 \rightarrow i-2, i+3 \rightarrow i-3, i+4 \rightarrow i-4$ y h$\rightarrow$ -h.
A continuación, procedemos a programar este algoritmo.
function deriv2num_h5(f,x)
f_2prima = zeros(length(f)) #vector de ceros del mismo tamaño que f, lo vamos a llenar
h = x[2] - x[1] #Los puntos siempre van a estar separados la distancia que haya entre el 1er y 2ndo elemento de x porque todos son equidistantes. (Malla uniforme)
doceh2 = 12*(h^2)
for i in 3:(length(f)-2) #el for corre sobre este intervalo ya que no estamos considerando los puntos extremos: 2 al inicio y 2 al final
f_2prima[i] = ( 16*(f[i+1]+ f[i-1]) - f[i+2] - f[i-2] - 30f[i] )/(doceh2)
end
#Calculo de los puntos extremos
#definimos los coeficientes
c0,c1,c2,c3,c4 = 35/12,-26/3,19/2,-14/3,11/12
#calculo del primer punto
f_2prima[1] = ( c0*f[1] + c1*f[2] + c2*f[3] + c3*f[4] + c4*f[5] )/h^2
# calculo del segundo punto
f_2prima[2] = ( c0*f[2] + c1*f[3] + c2*f[4] + c3*f[5] + c4*f[6] )/h^2
# calculo del ultimo punto o el punto N
f_2prima[end] = ( c0*f[end] + c1*f[end-1] + c2*f[end-2] + c3*f[end-3] + c4*f[end-4] )/h^2
# calculo del penultimo punto o el N-1
f_2prima[end-1] = ( c0*f[end-1] + c1*f[end-2] + c2*f[end-3] + c3*f[end-4] + c4*f[end-5] )/h^2
return f_2prima
end
A continuación probaremos nuestra función.
x = collect(range(0,stop=2π,length=50))
y = cos.(x)
y_2deriv = -cos.(x)
y_num_2deriv = deriv2num_h5(y,x)
plot(x,y,label="f(x)")
plot!(x,y_2deriv,label="f''(x)")
scatter!(x,y_num_2deriv,label="f''(x) numérica",m=:cross)
En este caso, vamos a obtener la derivada de una función en el punto $x_0$ pero con una precisión en la que se desprecian términos en $h^3$, i.e. tenemos que realizar operaciones sólo con un punto a la izquierda y uno a la derecha para una expresión simétrica (orden 2). Sin embargo, como se trata de una malla no uniforme, los puntos no estan a la misma distancia h como en el inciso anterior. Suponemos un punto a la izquierda de x_0 a distancia $h_1$ y uno a la derecha a distancia $h_2$.
Por lo tanto, partimos de las siguientes dos ecuaciones:
$$f(x_0-h_1)=f(x_0)-f'(x_0)h_1 + \frac{1}{2}f''(x_0)h_1^2 + \theta(h^3)\quad...(7)$$$$f(x_0+h_2)=f(x_0)+f'(x_0)h_2 + \frac{1}{2}f''(x_0)h_2^2 + \theta(h^3)\quad...(8)$$Así, podemos hacer $h_1^2$Eq(8) y $h_2^2$Eq(7) y obtenemos:
$$ h_2^2f(x_0-h_1)=h_2^2f(x_0)-f'(x_0)h_1h_2^2 + \frac{1}{2}f''(x_0)h_1^2h_2^2 + \theta(h^3)\quad...(9) $$$$ h_1^2f(x_0+h_2)=h_1^2f(x_0)+f'(x_0)h_1^2h_2 + \frac{1}{2}f''(x_0)h_2^2h_1^2 + \theta(h^3)\quad...(10) $$Restamos las ecuaciones (10)-(9):
$$ h_1^2f(x_0+h_2)-h_2^2f(x_0-h_1)=f(x_0)(h_1^2-h_2^2)+ f'(x_0)[h_1^2h_2+h_1h_2^2] $$Por lo que, la expresión para la derivada numérica de orden 2 en una malla no uniforme es:
$$ f'(x_0)=\frac{h_1^2f(x_0+h_2)-h_2^2f(x_0-h_1)-f(x_0)(h_1^2-h_2^2)}{h_1^2h_2+h_1h_2^2} $$Y recordando que supusimos que el punto a la izquierda de $x_0$ era el que estaba a distancia $h_1$ y el punto a la derecha a distancia $h_2$:
$$ \Rightarrow f'_i=\frac{h_1^2f_{i+1}-h_2^2f_{i-1}-f_i(h_1^2-h_2^2)}{h_1^2h_2+h_1h_2^2} $$Ahora, observamos que la expresión anterior no puede ser aplicada en los puntos extremos. Para calcular la derivada en los extremos necesitamos las expresiones para $x_0+h_1$ y $x_0+h_2$. Estas expresiones podemos escribirlas como el siguiente sistema de ecuaciones:
$$\begin{bmatrix} h_1 & \frac{1}{2}h_1^2 \\ h_2 & \frac{1}{2}h_2^2 \\ \end{bmatrix} \begin{bmatrix} f'(x_0) \\ f''(x_0) \end{bmatrix} = \begin{bmatrix} f(x_0+h_1)-f(x_0) \\ f(x_0+h_2)-f(x_0) \end{bmatrix} $$$\therefore$ La solución para $f(x_0)'$ es la solución del primer nivel y es la siguiente:
$$ f'(x_0)=\frac{h_2^2 \big (f(x_0+h_1)-f(x_0) \big )-h_1^2 \big (f(x_0 + h_2)-f(x_0) \big )}{h_1h_2^2-h_1^2h_2} $$$$ \Rightarrow f'(x_0)=\frac{h_2^2 f(x_0+h_1)- h_2^2f(x_0)-h_1^2 f(x_0 + h_2)+ h_1^2 f(x_0)}{h_1h_2^2-h_1^2h_2} $$En donde, para el primer punto se evalúa en $i=1$ y se tiene que $h_1$ corresponde a la distancia entre el primer y segundo punto, y $h_2$ a la distancia entre el primer y el 3er punto. ($i+1$ y $i+2$)
Por lo que, la expresión para el primer punto es la siguiente:
$$ f'_1=\frac{h_2^2 f_2- h_2^2f_1-h_1^2 f_3+ h_1^2 f_1}{h_1h_2^2-h_1^2h_2}$$Y para el último punto, se evalúa en $i=end$ y se tiene que $h_1$ corresponde a la distancia entre el último y el penúltimo punto, y $h_2$ a la distancia entre el último y el antepenúltimo punto. Por lo que, hay que realizar el siguiente cambio: $i+1 \rightarrow i-1, i+2 \rightarrow i-2$, y considerar un cambio de signo a todo.
Por lo tanto, la expresión para el último punto es la siguiente:
$$ f'_{end}=- \bigg ( \frac{h_2^2 f_{end-1}- h_2^2f_{end}-h_1^2 f_{end-2}+ h_1^2 f_{end}}{h_1h_2^2-h_1^2h_2} \bigg )$$A continuación, procedemos a programar este algoritmo.
function derivnum_h3(f,x)
f_prima = zeros(length(f)) #vector de ceros del mismo tamaño que f, lo vamos a llenar
for i in 2:(length(x)-1) #el for corre sobre este intervalo ya que no estamos considerando los puntos extremos
h1=x[i]-x[i-1] #h1 es la distancia entre el punto i y el punto a la izquierda
h2=x[i+1]-x[i] #h2 es la distancia entre el punto i y el punto a la derecha
f_prima[i]= (((h1^2)*f[i+1]-(h2^2)*f[i-1])-(f[i]*((h1^2)-(h2^2))))/(((h1^2)*h2)+(h1*(h2^2)))
end
#Calculo de los extremos
#calculo del primer punto
hder1=x[2]-x[1] #Para el primer punto, los dos puntos que va a consideran van a ser los dos a la derecha
hder2=x[3]-x[1]
f_prima[1]=(((hder2^2)*f[2])-((hder2^2)*f[1])-((hder1^2)*f[3])+((hder1^2)*f[1]))/((hder1*(hder2^2))-((hder1^2)*hder2))
#calculo del ultimo punto
hizq1=x[end]-x[end-1] ##Para el ultimo punto, los dos puntos que va a consideran van a ser los dos a la izquierda
hizq2=x[end]-x[end-2]
f_prima[end]=-(((hizq2^2)*f[end-1])-((hizq2^2)*f[end])-((hizq1^2)*f[end-2])+((hizq1^2)*f[end]))/((hizq1*(hizq2^2))-((hizq1^2)*hizq2))
return f_prima
end
x = collect(range(0,stop=2π,length=50))
y = cos.(x)
y_deriv = -sin.(x)
y_num_deriv = derivnum_h3(y,x)
plot(x,y,label="f(x)")
plot!(x,y_deriv,label="f'(x)")
scatter!(x,y_num_deriv,label="f'(x) numérica",m=:cross)
Vamos a obtener la derivada de una función en el punto $x_0$ con una precisión en la que se desprecian términos en $h^3$, i.e. tenemos que realizar operaciones sólo con un punto a la izquierda y uno a la derecha para una expresión simétrica (orden 2). Al igual que en el inciso anterior, se trata de una malla no uniforme, por lo que, suponemos un punto a la izquierda a distancia $h_1$ y uno a la derecha a distancia $h_2$.
Por lo tanto, podemos partir de las mismas dos ecuaciones que en b):
$$ f(x_0-h_1)=f(x_0)-f'(x_0)h_1 + \frac{1}{2}f''(x_0)h_1^2 + \theta(h^3)\quad...(7) $$$$ f(x_0+h_2)=f(x_0)+f'(x_0)h_2 + \frac{1}{2}f''(x_0)h_2^2 + \theta(h^3)\quad...(8) $$Asi, podemos hacer $h_1$Eq(8) y $h_2$Eq(7), y obtenemos:
$$ h_2f(x_0-h_1)=h_2f(x_0)-f'(x_0)h_1h_2 + \frac{1}{2}f''(x_0)h_1^2h_2 + \theta(h^3)\quad...(11) $$$$ h_1f(x_0+h_2)=h_1f(x_0)+f'(x_0)h_1h_2 + \frac{1}{2}f''(x_0)h_1h_2^2 + \theta(h^3)\quad...(12) $$Sumamos las ecuaciones (11) y (12):
$$ h_2f(x_0-h_1)+h_1f(x_0+h_2)=f(x_0)(h_1+h_2)+ \frac{1}{2}f''(x_0)[h_1h_2^2+h_1^2h_2] $$Por lo que, la expresión para la segunda derivada numérica de orden 2 en una malla no uniforme es:
$$ f''(x_0)= 2 \bigg (\frac{h_2f(x_0-h_1)+h_1f(x_0+h_2)-f(x_0)(h_1+h_2)}{[h_1h_2^2+h_1^2h_2]} \bigg) $$Y recordando que supusimos que el punto a la izquierda de $x_0$ era el que estaba a distancia $h_1$ y el punto a la derecha a distancia $h_2$:
$$ \Rightarrow f''_i = 2 \bigg (\frac{h_2f_{i-1}+h_1f_{i+1}-f_i(h_1+h_2)}{[h_1h_2^2+h_1^2h_2]} \bigg ) $$Ahora, observamos que la expresión anterior no puede ser aplicada en los puntos extremos. Para calcular la derivada en los extremos, como estamos en orden 2 al igual que en el inciso b), tenemos el mismo sistema de ecuaciones:
$$\begin{bmatrix} h_1 & \frac{1}{2}h_1^2 \\ h_2 & \frac{1}{2}h_2^2 \\ \end{bmatrix} \begin{bmatrix} f'(x_0) \\ f''(x_0) \end{bmatrix} = \begin{bmatrix} f(x_0+h_1)-f(x_0) \\ f(x_0+h_2)-f(x_0) \end{bmatrix} $$$\therefore$ La solución para $f(x_0)''$ es la solución del segundo nivel y es la siguiente:
$$f''(x_0)=2 \bigg (\frac{h_1 \big (f(x_0+ h_2)-f(x_0) \big )-h_2 \big (f(x_0+h_1)-f(x_0) \big )}{h_1h_2^2-h_1^2h_2} \bigg )$$$$\Rightarrow f''(x_0)=2\bigg ( \frac{h_1f(x_0 + h_2)-h_2f(x_0+h_1) +f(x_0)(h_2-h_1)}{h_1h_2^2-h_1^2h_2} \bigg )$$En donde, para el primer punto se evalúa en $i=1$ y se tiene que $h_1$ corresponde a la distancia entre el primer y segundo punto, y $h_2$ a la distancia entre el primer y el 3er punto. ($i+1$ y $i+2$)
Por lo que, la expresión para el primer punto es la siguiente:
$$ f''_1=2\bigg ( \frac{h_1f_3-h_2f_2 +f_1(h_2-h_1)}{h_1h_2^2-h_1^2h_2} \bigg )$$Y para el último punto, se evalúa en $i=end$ y se tiene que $h_1$ corresponde a la distancia entre el último y el penúltimo punto, y $h_2$ a la distancia entre el último y el antepenúltimo punto. Por lo que, hay que realizar el siguiente cambio: $i+1 \rightarrow i-1, i+2 \rightarrow i-2$.
Por lo tanto, la expresión para el último punto es la siguiente:
$$ f''_{end}=2\bigg ( \frac{h_1f_{end-2}-h_2f_{end-1} +f_{end}(h_2-h_1)}{h_1h_2^2-h_1^2h_2} \bigg )$$A continuación procedemos a programar este algoritmo.
function deriv2num_h3(f,x)
f_2prima = zeros(length(f))
for i in 2:(length(f)-1)
h1=x[i]-x[i-1] #h1 es la distancia entre el punto i y el punto a la izquierda
h2=x[i+1]-x[i] #h2 es la distancia entre el punto i y el punto a la derecha
f_2prima[i] = 2*( ((h2*f[i-1])+(h1*f[i+1])-(f[i]*(h1+h2)))/((h1*(h2^2))+((h1^2)*h2)) )
end
#Calculo de los extremos
#calculo del primer punto
hder1=x[2]-x[1] #Para el primer punto, los dos puntos que va a consideran van a ser los dos a la derecha
hder2=x[3]-x[1]
f_2prima[1] = 2*( ((hder1*f[3])-(hder2*f[2])+(f[1]*(hder2-hder1)))/((hder1*(hder2^2))-((hder1^2)*hder2)) )
#calculo del ultimo punto
hizq1=x[end]-x[end-1] ##Para el ultimo punto, los dos puntos que va a consideran van a ser los dos a la izquierda
hizq2=x[end]-x[end-2]
f_2prima[end] = 2*( ((hizq1*f[end-2])-(hizq2*f[end-1])+(f[end]*(hizq2-hizq1)))/((hizq1*(hizq2^2))-((hizq1^2)*hizq2)) )
return f_2prima
end
x = collect(range(0,stop=2π,length=50))
y = cos.(x)
y_2deriv = -cos.(x)
y_num_2deriv = deriv2num_h3(y,x)
plot(x,y,label="f(x)")
plot!(x,y_2deriv,label="f''(x)")
scatter!(x,y_num_2deriv,label="f''(x) numérica",m=:cross)
Importamos primero la biblioteca que utilizaremos para medir tiempos.
using BenchmarkTools
Ahora, las funciones que programamos en el problema anterior, calculan la primera o segunda derivada de una función de manera numérica, por lo que, podemos escoger cualquier función para probar el tiempo que se tarda. Escogemos sin(x).
x = collect(1:1:500) #Definimos una malla con 500 puntos
y = sin.(x) #Funcion a la cual le vamos a calcular la derivada para cada punto
println("Número de puntos: ",length(x))
@benchmark deriv2num_h5(y,x)
$\therefore$ El tiempo promedio es 4.684 $\mu$s.
#Construimos una malla no uniforme de 500 puntos
x_a = collect(1:2:500) #De 1 a 500 en pasos de 2
x_b = collect(501:1:750) #De 501 a 750 en pasos de 1
x = vcat(x_a,x_b) #malla compuesta
y = sin.(x) #Funcion a la cual le vamos a calcular la derivada para cada punto
println("Número de puntos: ",length(x))
@benchmark derivnum_h3(y,x)
$\therefore$ El tiempo promedio es 5.698 $\mu$s.
#Construimos una malla no uniforme de 500 puntos
x_a = collect(1:2:500) #De 1 a 500 en pasos de 2
x_b = collect(501:1:750) #De 501 a 750 en pasos de 1
x = vcat(x_a,x_b) #malla compuesta
y = sin.(x) #Funcion a la cual le vamos a calcular la derivada para cada punto
println("Número de puntos: ",length(x))
@benchmark deriv2num_h3(y,x)
$\therefore$ El tiempo promedio es 4.633 $\mu$s.
Para el caso de la segunda derivada tenemos las funciones definidas en a) y c), observamos que la de a) tiene 2 órdenes más de precisión y además es casi igual de rápida que la de c), por lo tanto, presenta dos ventajas sobre ésta ya que es mucho más eficiente, sin embargo, su algoritmo fue el más complicado de todos y sólo funciona para mallas uniformes. Luego, por su parte, el algoritmo de c) presenta la ventaja de que es mucho más rápido de programar debido a que el sistema de ecuaciones que se tuvo que resolver para encontrar la expresión era mucho más sencillo. Asimismo, incluye mallas no uniformes, por lo que se trata de un caso más general que el de a).
Por otro lado, en el caso de la primera derivada tenemos la función definida en b), en este caso, aunque es el algoritmo más lento de los 3, tampoco podemos decir que se tarde mucho más que los otros, además es rápido de programar y considera mallas no uniformes al igual que el de c).
x1=collect(range(-5,stop=5,length=600)) #intervalo de la x: de -5 a 5
f=exp.(-(x1.^2)).*sin.(10*x1) #definicion de f
df= (2*exp.(-(x1.^2))).*((5*cos.(10*x1)) - (x1.*sin.(10*x1))) #derivada analítica de f
d2f= (exp.(-(x1.^2))).*(((2*((2*(x1.^2)).-51)).*(sin.(10*x1))).-((40*x1).*cos.(10*x1))) #segunda derivada analítica de f
plot(x1,f,label="f(x)") #Grafica de f
plot!(x1,df,label="f'(x)") #Grafica de f'
plot!(x1,d2f,label="f''(x)") #Grafica de f''
x2=collect(range(0.5,stop=5,length=500)) #intervalo de las x: de 0.5 a 5
g = -2*((x2.^6).^(-1)) + 1*((x2.^(12)).^(-1)) #definicion de g
dg = (12*(x2.^(6).-1)).*(x2.^(-13)) #derivada analítica de g
d2g =(-84*(x2.^6).+156).*(x2.^(-14)) #segunda derivada analítica de g
plot(x2,g,label="g(x)") #Grafica de g
plot!(x2,dg,label="g'(x)") #Grafica de g'
plot!(x2,d2g,label="g''(x)") #Grafica de g''
Variamos el tamaño del paso o la longitud del collect hasta encontrar el adecuado. (Solo pondremos el código de la malla adecuada)
a) Segunda derivada de orden 4 para malla uniforme.
$f(x)=e^{-x^2} \sin(10x)$:
#Funcion f(x)
x1 = collect(range(-5,stop=5,length=1949)) #intervalo de la x: de -5 a 5
f=exp.(-(x1.^2)).*sin.(10*x1) #definicion de f
d2f= (exp.(-(x1.^2))).*(((2*((2*(x1.^2)).-51)).*(sin.(10*x1))).-((40*x1).*cos.(10*x1))) #segunda derivada analítica de f
y2_h5 = deriv2num_h5(f,x1) #segunda derivada numérica de orden 4 de f(x)
println("Error con h5 = " , maximum(abs.(d2f-y2_h5)) ) #Imprime el error maximo entre y2_h5 y d2f
$\therefore$ La malla adecuada para que el error sea menor a $10^{-5}$ es: x1 = collect(range(-5,stop=5,length=1949))
$g(x)=-\frac{2}{x^6} + \frac{1}{x^{12}}$:
#Funcion g(x)
x2=collect(range(0.5,stop=5,length=90000)) #intervalo de las x: de 0.5 a 5
g = -2*((x2.^6).^(-1)) + 1*((x2.^(12)).^(-1)) #definicion de g
d2g =(-84*(x2.^6).+156).*(x2.^(-14)) #segunda derivada analítica de g
y2_h5 = deriv2num_h5(g,x2) #segunda derivada numérica de orden 4 de g(x)
println("Error con h5 = " , maximum(abs.(d2g-y2_h5)) ) #Imprime el error maximo entre y2_h5 y d2g
En este caso, observamos que el mínimo error que pudimos conseguir con nuestro algoritmo fue del orden de $10^{-2}$, ésto debido a que la función g(x) a diferencia de la f(x) tiene un comportamiento extraño justo en el intervalo que estamos considerando, como podemos apreciar en la gráfica que está al incio de este problema, ya que, específicamente cerca del 0.5 su segunda derivada diverge. Y como, en este algoritmo se considera la serie de Taylor despreciando términos de $h^5$, podemos concluir que esa precisión no es suficiente para este caso.
b) Derivada de orden 2 para malla no uniforme
$f(x)=e^{-x^2} \sin(10x)$:
#Funcion f(x)
#Construimos una malla no uniforme
x_a=collect(-5:0.00002:-4)
x_b=collect(-4.1:0.0001:5)
x1=vcat(x_a,x_b)
f=exp.(-(x1.^2)).*sin.(10*x1) #definicion de f
df= (2*exp.(-(x1.^2))).*((5*cos.(10*x1)) - (x1.*sin.(10*x1))) #derivada analítica de f
y_h3 = derivnum_h3(f,x1) #derivada numérica de orden 2 de f(x)
println("Error con h3 = " , maximum(abs.(df-y_h3)) ) #Imprime el error maximo entre y_h3 y df
$\therefore$ La malla adecuada para que el error sea menor a $10^{-5}$ es: `
x_a=collect(-5:0.00002:-4)
x_b=collect(-4.1:0.0001:5)
x1=vcat(x_a,x_b)
$g(x)=-\frac{2}{x^6} + \frac{1}{x^{12}}$:
#Funcion g(x)
#Construimos una malla no uniforme
x_a= collect(0.5:0.000001:4)
x_b= collect(3.9:0.0000001:5)
x2= vcat(x_a,x_b)
g = -2*((x2.^6).^(-1)) + 1*((x2.^(12)).^(-1)) #definicion de g
dg = (12*(x2.^(6).-1)).*(x2.^(-13)) #derivada analítica de g
y_h3 = derivnum_h3(g,x2) #derivada numérica de orden 2 de g(x)
println("Error con h3 = " , maximum(abs.(dg-y_h3)) ) #Imprime el error maximo entre y_h3 y dg
$\therefore$ La malla adecuada para que el error sea del orden de $10^{-5}$ es:
x_a= collect(0.5:0.000001:4)
x_b= collect(3.9:0.0000001:5)
x2= vcat(x_a,x_b)
c) Segunda derivada de orden 2 para malla no uniforme
$f(x)=e^{-x^2} \sin(10x)$:
#Funcion f(x)
#Construimos una malla no uniforme
x_a= collect(-5:0.00004:-4)
x_b= collect(-4.1:0.00005:5)
x1= vcat(x_a,x_b)
f=exp.(-(x1.^2)).*sin.(10*x1) #definicion de f
d2f= (exp.(-(x1.^2))).*(((2*((2*(x1.^2)).-51)).*(sin.(10*x1))).-((40*x1).*cos.(10*x1))) #segunda derivada analítica de f
y2_h3 = deriv2num_h3(f,x1) #segunda derivada numérica de orden 2 de f(x)
println("Error con h3 = " , maximum(abs.(d2f-y2_h3)) ) #Imprime el error maximo entre y2_h3 y d2f
$\therefore$ La malla adecuada para que el error sea menor a $10^{-5}$ es: `
x_a= collect(-5:0.00004:-4)
x_b= collect(-4.1:0.00005:5)
x1= vcat(x_a,x_b)
$g(x)=-\frac{2}{x^6} + \frac{1}{x^{12}}$:
# Funcion g(x)
#Construimos una malla no uniforme
x_a=collect(range(0.5,stop=1.1,length=1000000))
x_b=collect(range(1,stop=5,length=1000))
x2 = vcat(x_a,x_b)
g = -2*((x2.^6).^(-1)) + 1*((x2.^(12)).^(-1)) #definicion de g
d2g =(-84*(x2.^6).+156).*(x2.^(-14)) #segunda derivada analítica de g
y2_h3 = deriv2num_h3(g,x2) #segunda derivada numérica de orden 2 de g(x)
println("Error con h3 = " , maximum(abs.(d2g-y2_h3)) ) #Imprime el error maximo entre y2_h3 y d2g
En este caso, observamos que no pudimos conseguir un error de $10^{-5}$, ésto debido a las mismas razones planteadas en el inciso a) ya que también estamos obteniendo la segunda derivada de g(x). Sin embargo, a diferencia de a), el error en este caso es del orden de $10^2$, ya que se trata de un algoritmo de orden 2, es decir, estamos despreciando términos del orden de $h^3$, y por tanto, su precisión se reduce todavía más en comparación con el algoritmo de a).
Las mallas no siempre tienen el mismo espacio entre sus puntos, recordemos que en el inciso a) el algoritmo fue programado pensando en mallas uniformes por lo que en dicho caso, si se utilizó una malla uniforme. Sin embargo, para los incisos b) y c), los algoritmos fueron pensados para mallas no uniformes, y justamente, en esos casos, encontramos mallas no uniformes. Es importante mencionar que a pesar de que en la mayoría de los casos, si encontramos la malla adecuada para que el error fuera de $10^{-5}$, tuvimos que utilizar mallas con una cantidad de puntos muy grandes, lo que revela que la precisión de nuestros algoritmos no es la más óptima.
Recordamos el algortimo simétrico visto en clase para la derivada numérica de orden 2 para una malla uniforme:
function deriv_h3(f,x)
f_prima = zeros(length(f))
h = x[2] - x[1]
dosh = 2*h
for i in 2:(length(f)-1)
f_prima[i] = ( f[i+1]-f[i-1] )/(dosh)
end
#Calculo de los extremos
# calcular el primer punto
f_prima[1] = ( -0.5*f[3] + 2*f[2] -1.5*f[1] )/h
# calcular el ultimo punto
f_prima[end] = -( -0.5*f[end-2] + 2*f[end-1] -1.5*f[end] )/h
return f_prima
end
Y el algoritmo simétrico para la derivada numérica de orden 4 para una malla uniforme también visto en clase:
function deriv_h5(f,x)
f_prima = zeros(length(f))
h = x[2] - x[1]
doceh = 12*h
for i in 3:(length(f)-2)
f_prima[i] = ( -f[i+2]+ 8*(f[i+1]-f[i-1]) +f[i-2] )/(doceh)
end
#Calculo de los extremos
c0,c1,c2,c3,c4 = -25/12,4,-3,4/3,-0.25
# calcular el primer punto
f_prima[1] = ( c0*f[1] + c1*f[2] + c2*f[3] + c3*f[4] + c4*f[5] )/h
# calcular el segundo punto
f_prima[2] = ( c0*f[2] + c1*f[3] + c2*f[4] + c3*f[5] + c4*f[6] )/h
# calcular el ultimo punto o el punto N
f_prima[end] = -( c0*f[end] + c1*f[end-1] + c2*f[end-2] + c3*f[end-3] + c4*f[end-4] )/h
# calcular el penultimo punto o el N-1
f_prima[end-1] = -( c0*f[end-1] + c1*f[end-2] + c2*f[end-3] + c3*f[end-4] + c4*f[end-5] )/h
return f_prima
end
Tiempos de ejecución para calcular la derivada utillizando el algoritmo de 2ndo orden:
x = collect(-2π:1:2π) #Intervalo de las x: -2π a 2π (Malla uniforme)
sinc = (sin.(3*x)).*(3*x.^(-1)) #definicion de la funcion sinc(x) (Funcion a la cual le vamos a calcular la derivada para cada punto)
println("Número de puntos: ",length(x))
@benchmark deriv_h3(sinc,x)
Tiempos de ejecuación para calcular la derivada utilizando el algoritmo de 4to orden:
x = collect(-2π:1:2π) #Intervalo de las x: -2π a 2π (Malla uniforme)
sinc = (sin.(3*x)).*(3*x.^(-1)) #definicion de la funcion sinc(x) (Funcion a la cual le vamos a calcular la derivada para cada punto)
println("Número de puntos: ",length(x))
@benchmark deriv_h5(sinc,x)
Variamos el tamaño del paso dentro del collect hasta encontrar el adecuado para que el error sea menor a $10^{-5}$. (Solo pondremos el código del tamaño del paso adecuado)
Algoritmo de orden 2
x=collect(-2π:0.001:2π) #intervalo de x/ el paso en este caso es 0.001
sinc=(sin.(3*x)).*(3*x.^(-1)) #definicion de sinc(x)
Dsinc= (((3*x).*(cos.(3*x))).-(sin.(3*x))).*(3*x.^(-2)); #derivada analítica de sinc(x)
y_h3 = deriv_h3(sinc,x) #derivada numérica de orden 2 de sinc(x)
println("Error con h3 =" , maximum(abs.(Dsinc-y_h3)) ) #Imprime el error maximo entre y_h3 y Dsinc
plot(x,Dsinc,label="sinc''(x) analítica")
plot!(x,y_h3,label="sinc''(x) numérica orden 2",m=:cross)
$\therefore$ El espacio entre los puntos de la malla para que el error máximo sea menor a $10^{-5}$ para el algoritmo de 2ndo orden es: 0.001
Algoritmo de orden 4
x=collect(-2π:0.02:2π) #intervalo de x/ el paso en este caso es 0.02
sinc=(sin.(3*x)).*(3*x.^(-1)) #definicion de sinc(x)
Dsinc= (((3*x).*(cos.(3*x))).-(sin.(3*x))).*(3*x.^(-2)); #derivada analítica de sinc(x)
y_h5 = deriv_h5(sinc,x) #derivada numérica de orden 4 de sinc(x)
println("Error con h5 =" , maximum(abs.(Dsinc-y_h5)) ) #Imprime el error máximo entre y_h5 y Dsinc
plot(x,Dsinc,label="sinc''(x) analítica")
plot!(x,y_h5,label="sinc''(x) numérica orden 4",m=:cross)
$\therefore$ El espacio entre los puntos de la malla para que el error máximo sea menor a $10^{-5}$ para el algoritmo de 4to orden es: 0.02
Vamos a ver primero cuantos puntos tiene una malla x con un espaciamiento de 0.25:
x=collect(-2π:0.25:2π)
length(x)
Ahora si, para cada algoritmo vemos igual cuantos puntos tiene la malla con el espaciamiento encontrado en el inciso anteior, y después procedemos a graficar.
Algoritmo de orden 2
x_1=collect(-2π:0.001:2π) #malla encontrada en el inciso anterior
length(x_1)
Por lo tanto, partimos de 51 puntos hasta 12567.
num_puntos = collect(51:596:12567) #Tomamos 596 como paso porque 12567-51=12516 y 12516/596=21
#i.e. estamos considerando 21+1 (porque Julia si considera el primero) puntos para graficar.
tiempos_promedio = zeros(length(num_puntos))
for n_i in 1:length(num_puntos)
x = collect(range(51,step=596,length=num_puntos[n_i]))
sinc=(sin.(3*x)).*(3*x.^(-1))
y_h3=deriv_h3(sinc,x)
t = @benchmark deriv_h3(sinc,x)
tiempos_promedio[n_i] = mean(t.times)
end
plot(num_puntos,tiempos_promedio,
marker=:dot,legend=false,
title="Tiempos promedio para evaluar ",
xlabel="Número de puntos",ylabel="Tiempo en [ns]")
Algoritmo de orden 4
x_2=collect(-2π:0.02:2π) #malla encontrada en el inciso anterior
length(x_2)
Por lo tanto, partimos de 51 puntos hasta 629.
num_puntos = collect(51:34:629) #Tomamos 34 como paso porque 629-51=578 y 578/34=17
#i.e. estamos considerando 17+1 (porque Julia si considera el primero) puntos para graficar.
tiempos_promedio = zeros(length(num_puntos))
for n_i in 1:length(num_puntos)
x = collect(range(51,step=34,length=num_puntos[n_i]))
sinc=(sin.(3*x)).*(3*x.^(-1))
y_h5=deriv_h5(sinc,x)
t = @benchmark deriv_h5(sinc,x)
tiempos_promedio[n_i] = mean(t.times)
end
plot(num_puntos,tiempos_promedio,
marker=:dot,legend=false,
title="Tiempos promedio para evaluar ",
xlabel="Número de puntos",ylabel="Tiempo en [ns]")
Pudimos ver que el tiempo de ejecución en ambos algoritmos es casi el mismo, sin embargo, para el algoritmo de orden 4 necesitamos un orden de magnitud menor de espaciamiento en la malla, y por lo mismo mucho menos puntos (419 a 12567) para alcanzar un error menor a $10^{-5}$, por lo que el algoritmo de orden 4 nos ofrece mucho más precisión y se tarda la misma cantidad de tiempo en ejecutarse, razón por la que podemos decir que es mucho mejor que el algoritmo de orden 2. Por otra parte, tras realizar las gráficas del inciso anterior, observamos que el tiempo de ejecución va aumentando linealmente (aproximadamente) conforme aumenta el número de puntos en la malla, lo cual es lógico, ya que tiene sentido pensar que entre más puntos tenga, más operaciones realizará y por tanto se tardará más tiempo.